#include <iostream>
#include <boost\shared_ptr.hpp>
#include <boost\progress.hpp>
#include <boost\random\linear_congruential.hpp>
#include <boost\random\uniform_01.hpp>

#include "VarianceGammaDistribution.h"
#include "VGQAsympExp.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "Polynomial.h"
#include "VGFAsympExp.h"
#include "VGQLeftSeries.h"
#include "VGQRightSeries.h"
#include "VGQApprox.h"
#include "NormalQApprox.h";

#include "VGQSeries.h"
#include <iomanip>

int main() {
  
	using namespace std;
	
	double mu = 0.0002585;
	double lambda = 2.262443438914027;
	double beta = -2.3421746548782383;
	double alpha = 264.9366251164119;

	/*double mu = 20;
	double lambda = 6.556247;
	double beta = 15.335646548782383;
	double alpha = 25.96471164119;*/

	double accuracy = 1e-10;

	try{
		cout << "setup time: " << endl;
		shared_ptr<VGQApprox> pvgQapprox;
		{
			boost::progress_timer t; 
			pvgQapprox = make_shared<VGQApprox>(mu, lambda, alpha, beta, accuracy);
		}
		string outputDirectory = "C:\\WorkArea\\PhD\\TestProjects\\VC++\\QuantileMechanics\\ProfileResults\\VG"; 
		
		cout<< "Profiling VG quantile. Time elapsed: " << endl;

		{
			boost::progress_timer t; 
			pvgQapprox->profile(outputDirectory);
		}

		cout << "See output directory " << outputDirectory << " for results." << endl;
	
		boost::minstd_rand intgen;
		boost::uniform_01<boost::minstd_rand> gen(intgen);

		cout<< "1 million U[0,1] variates generated in: " << endl;
		const size_t n = 1e6;
		vector<double> uVariates(n);
		{
			boost::progress_timer t; 
			for(size_t i=0; i<n; ++i)
			{
				uVariates[i] = gen();
			}
		}
		
		cout<< "1 million VG quantile evaluations took: " << endl;

		{
			boost::progress_timer t; 
			for(int k = 0; k < n; k++)
			{
				(*pvgQapprox)(uVariates[k]);
			}
		
		}

		cout<< "1 million standard normal quantile evaluations took: " << endl;
		NormalQApprox normalQuantile(0.23563, 0.0535);
		{
			boost::progress_timer t; 
			for(int k = 0; k < n; k++)
			{
				normalQuantile(uVariates[k]);
			}
		}
	}catch(const std::exception& e){
		cout << "error: " << e.what() << endl;
	}

	int i;
    std::cin >> i;

	return 0;
}
